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CHAPTER  I 


INTRODUCTION 

1.1  Exponential  Analysis  as  a  Special  Case 
of  Systems  Identification 

Exponential  analysis  attempts  to  characterize  a  waveform  with  a 
sum  of  complex  exponentials,  that  is,  a  s* m  of  damped  sinusoidal  com¬ 
ponents-  Consider  the  class  of  linear  processes  whose  impulse  re¬ 
sponses  are  representable  as  a  sum  of  exponential  components.  If  the 
impulse  response,  possibly  noise  contamined,  is  given  for  a  process  in 
this  class,  the  transfer  function  of  that  process  can  be  estimated  by 
applying  a  body  of  theory  known  as  systems  identification  [I, 2, 3, 4]. 
The  impulse  response  can  be  expressed  as  the  inverse  transform  (either 
the  inverse  Laplace  transform  or  the  inverse  z-transform)  of  the 
partial-fraction  expansion  of  this  estimated  transfer  function.  If 
the  waveform  to  be  analyzed  is  assumed  to  be  the  impulse  response 
of  a  process  of  this  class,  then  the  systems  identification  technique 
plus  the  process  of  partial-fraction  expansion  can  be  viewed  as  an 
exponential  analysis  method.  Hence,  exponential  analysis  methods  can 
be  equated  to  systems  identification  methods  for  the  case  of  impulse 
input  to  the  process.  The  poles  of  the  transfer  function  are  the 
damped  resonances  that  characterize  the  waveform.  The  imaginary 
parts  of  the  poles  are  the  angular  frequencies  of  the  sinusoidal  com¬ 
ponents  and  the  real  parts  are  the  corresponding  damping  constants. 
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1.2  Motivation  for  Exponential  Analysis  and 
a  Survey  of  the  Literature 

Systems  identification  theory  and  exponential  analysis  find  appli¬ 
cations  in  such  diverse  fields  as  industrial  controls,  economic  model¬ 
ing  and  in  the  analysis  of  biological  systems.  Recently,  these  iden¬ 
tification  methods  have  found  application  in  the  extraction  of  the 
singularity  expansion  method  (SEM)  description  of  a  transient  scatterer 
from  its  time  domain  response  as  was  first  suggested  by  Mittra  and  Van 
Blaricum  [5].  SOI  was  developed  by  Baum  fiS,  7]  from  the  insight  that 
the  transient  response  of  a  scatterer  resembles  a  sum  of  exponentially 
damped  sinusoids.  The  least-squares  Prony’s  method  was  proposed  by 
Van  Blaricum  and  Mittra  [8,  9]  as  a  means  of  obtaining  the  SEM  descrip¬ 
tion  from  the  transient  response  of  a  scatterer.  Pearson  and  Roberson 


[10]  have  since  developed  and  documented  a  method  of  obtaining  the  com¬ 
plete  SEM  description  of  a  scatterer  from  transient  response  data. 
Dudley  [11]  related  Prony’s  method  to  a  parametric  system  model  and 
proceeded  to  demonstrate  a  bias  in  the  estimates  of  the  system  poles 
inherent  in  least-squares  Prony's  method. 

The  parametric  model  employed  by  Dudley  is  a  modified  version  of 
the  generalized  ^odel  described  by  Evkhoff  [1,  2],  Astrom  and  Eyknoff 
[3],  and  on  pages  209-220  of  Eyknoff  [4].  The  origin  of  the  general¬ 
ized  model  can  be  traced  to  Kalman  in  1958  [12]  who  assumes  noise-free 
input  and  output  records  of  the  process  to  be  identified.  This  is  the 
assumption  from  which  the  generalized  model  derives  its  validity.  With 
noise,  this  model  is  no  longer  valid,  and  the  resulting  transfer  func¬ 
tion  estimate  is  slightly  erroneous. 
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Steiglitz  and  McBride  [13]  introduced  a  so-called  actual  model 
that  differs  from  the  Kalman  or  generalized  model  by  using  the  model 
output  in  place  of  the  noise  corrupted  process  output  for  the  feedback 
inherent  in  the  model.  The  validity  of  the  actual  model  does  not  break 
down  when  noise  is  present  in  the  process  output.  However,  the  estima¬ 
tion  of  the  parameters  of  the  actual  model  is  a  highly  nonlinear  prob¬ 
lem.  One  major  reason  for  the  generalized  model’s  popularity  is  lin¬ 
earity  in  the  parameters  allowing  one-shot  estimation.  In  contrast,  an 
iterative  estimation  procedure  is  required  by  the  actual  model. 

Another  related  method  is  the  pencil-of-f unctions  method  advocated 
by  Sarkar  [14,  15].  The  method  was  originally  proposed  by  Jain  and 
Gupta  [16]  and  elaborated  on  by  Jain  [17,  18,  19]-  In  [15]  Sarkar 
indicates  connections  between  Prony’s  method,  the  Wiener  filter,  and 
the  pencil-of-functions  method. 

One  important  property  of  the  actual  model  is  that  when  it  is 
used  for  exponential  analysis,  it  produces  a  "best  fit"  to  the  wave¬ 
form  under  analysis  in  the  mean-square  3ense,  that  is,  it  minimizes 
the  mean-square  error.  Some  other  methods  that  have  this  property  are 
found  in  references  [13,20,21,22,23]. 

1.3  The  Contribution  of  the  Present  Work 

The  original  intent  of  this  work  was  to  develop  a  noise  tolerant, 
efficient  method  for  exponential  analysis.  The  method  was  to  find 
direct  application  in  extraction  of  the  SEM  description  of  a  scatterer 
from  measured  surface  currents.  A  new  noise  tolerant  method  is  presen¬ 
ted  in  this  document.  Unfortunately,  the  method  is  laborious,  and 
hence,  only  partial  success  can  be  claimed  with  regard  to  the  original 


intent. 


Perhaps  the  real  contribution  of  the  present  work  is  the  conceptu¬ 
al  groundwork  which  is  a  prerequisite  for  further  improvements  on  ex¬ 
ponential  analysis  methods.  This  groundwork  is  laid  by  defining  a 
general  scheme  that  incorporates  most,  if  not  all,  exponential  analysis 
methods  and  by  defining  the  sources  of  difficulties  for  such  methods. 
This  groundwork  is  established  through  a  new  formulation  of  an  extended 
generalized  model  and  an  extended  actual  model  after  Steiglitz  and 
McBride  [13] .  The  extension  is  to  admit  any  set  of  linear  filter 
functions  to  form  generalized  filter  sections  in  the  respective  models. 
Emphasis  is  placed  on  concept  development  rather  than  strict  mathe¬ 
matical  rigor  which  would  stifle  such  development. 

In  Chapter  II  one  major  source  of  difficulty  with  exponential 
analysis  methods  is  defined.  This  difficulty  is  the  parameter  bias 
that  Dudley  [11]  brought  to  light.  An  attempt  is  made  to  relate  the 
source  of  this  so-called  bias  to  the  distinction  made  by  Steiglitz 
and  McBride  [13]  between  the  true  error  model  and  the  linear  regres¬ 
sion  model.  Also  presented  in  Chapter  II  is  the  general  identifica¬ 
tion  scheme  that  is  used  through  the  rest  of  this  work. 

Chapter  III  relates  Prony's  method  and  the  pencil-of-f unctions 
method  to  the  general  scheme  of  Chapter  II.  The  various  problems  asso¬ 
ciated  with  the  pencil-of-f unctions  method  are  discussed,  and  simple 
remedies  to  those  problems  are  proposed. 

Chapter  IV  introduces  a  new  method,  called  the  adaptive  method, 
which  is  highly  tolerant  of  the  presence  of  noise  in  the  waveform 
under  analysis. 

Chapter  V  presents  conclusions  and  indicates  directions  in  which 
future  research  may  be  fruitful. 


1 . 4  Rotational  Conventions 


The  following  conventions  of  notation  are  observed: 

1.  The  symbol  "s"  denotes  the  Laplace-transform  variable.  If  the 
transfer  function  of  a  process  is  a  function  of  s,  the  process  is 
assumed  to  be  a  continuous-time  process. 

2.  The  symbol  "z"  denotes  the  z-transform  variable.  If  the  transfer 
function  of  a  process  is  a  function  of  z,  the  process  is  assumed 
to  be  a  discrete- time,  sampled-data  process. 

3.  If  a  symbol  for  a  transfer  function  appears  alone,  without  the 
appropriate  variable  enclosed  in  parenthesis,  the  symbol  denotes 
a  "generalized"  transfer  function  whose  functional  dependence  is 
not  restricted  in  any  way.  To  illustrate,  a  process  can  be  given  a 
c' ntinuous-time  representation  in  which  case  the  transfer  function 
could  be  a  function  of  s,  or  the  process  can  be  given  a  discrete¬ 
time  representation,  then  the  transfer  function  is  a  function  of  z. 
If  only  the  essence  of  the  process  is  to  be  conveyed,  without  re¬ 
gard  to  its  particular  representation,  then  the  transfer  function 
is  written  without  indicating  its  functional  dependence. 

4.  An  asterisk  indicates  the  complex  conjugate  of  the  expression  pre¬ 
ceding  it. 

5.  The  symbol  "T"  raised  above  the  expression  it  follows  indicates  the 
transpose  of  the  matrix  or  vector  preceding  it. 

6.  A  prime  indicates  the  transpose  conjugate  of  the  matrix  or  vector 
preceding  it. 

Other  symbols  used  in  this  work  are  defined  when  they  are  introduced 

in  the  text. 
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CHAPTER  II 


ESTIMATION  OF  THE  PROCESS  TRANSFER  FUNCTION 


2.1  Introduction 

The  problem  under  consideration  is  the  characterization  of  a  sam¬ 
pled,  noise  contaminated  waveform,  y(k) ,  as  a  weighted  sum  of  complex 
exponentials  of  the  form 

n  . 

y_(k)  *  £  A  z.,  k  -  0,  M-l  (1) 
j-1  J  3 

where  z^  =  exp(s^.T)  and  T  is  the  time  duration  between  successive 
samples.  The  objective  is  to  choose  the  A  and  z^  that  best  character¬ 
ize  the  waveform.  Usually  the  best  characterization  is  assumed  to  be 
the  one  that  minimizes  the  mean-squared  error  between  the  waveform  and 
the  characterization. 

It  is  assumed  that  y(k)  is  the  noise  contaminated,  impulse  re¬ 
sponse  of  the  linear  process  with  a  nth  order  transfer  function  given  by 


blPl  + 


+  b  F 
n  n 


a  +  a  F  + 
o  11 


+  a  F 
n  n 


(2) 


The  F^  are  predefined  functions  of  the  transform  variable.  For  example, 
the  could  be  polynomials  in  the  transform  variable.  Later  the  F.^ 
are  related  to  filtering  operations,  and  for  this  reason,  they  are  re¬ 
ferred  to  as  "filter  transfer  functions."  With  appropriate  choices  of 
the  e^,  b^,  and  F^,  any  transfer  function  can  be  represented  in  this 
form.  The  reason  this  form  is  used  is  because  the  general  parametric 
model  introduced  in  the  next  section  also  has  a  transfer  function  of 
this  form.  The  F.^  are  the  transfer  functions  (which  are  not  specified 
explicitly  in  order  to  remain  completely  general)  of  the  filters  that 
make  up  the  model. 
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The  problem  can  be  restated  as  the  estimation  of  (2).  By  expand¬ 
ing  the  estimated  transfer  function  in  partial  fractions  and  performing 
an  inverse  transform,  the  desired  characterization  is  obtained.  The  z ^ 
or  can  be  interpreted  as  the  poles  of  (2)  and  the  A^  are  the  corres¬ 
ponding  residues  of  (2).  The  z^  and  s^  are  related  by  z^  *  exp(siT). 

2.2  Two  Parametric  Models 

Steiglitz  and  McBride  describe  two  parametric  models  in  terms  of 
which  they  interpret  their  identification  procedure.  A  straightfor¬ 
ward  extension  of  these  parametric  models  provides  a  framework  that  is 
broader  in  its  applicability.  They  are  specialized  herein  to  describe 
the  Prony  procedure  [24],  the  so-called  least  squares  Prony  procedure 
[9],  and  the  pencil-of-functions  method  of  Jain  [19].  The  models  are 
used  further  to  provide  a  new  iterative  identification  procedure  that 
appears  to  be  substantially  more  tolerant  to  noise-corruption  in  data 
than  existing  schemes.  The  first,  the  actual  or  true  model,  is  shown 
in  Figure  1  and  has  the  transfer  function 

S.F  +  •••  +  8  F 
u  —  11  n  n _ 

m  a  +  a.F .  +  •••  +  a  F  * 
oil  n  n 

The  second  model  is  called  the  approximate  model  and  is  shown  in 
Figure  2.  The  approximate  model  is  derived  from  the  actual  model  by 
approximating  the  model  output  with  the  noisy  process  output  as  the 
source  of  feedback  in  the  model.  It  should  be  noted  that  (3)  is  not 
the  transfer  function  of  the  approximate  model.  In  fact,  a  transfer 
function  that  serves  as  a  useful  approximation  to  the  process  transfer 
function  cannot  be  defined  for  the  approximate  model.  These  models 
reduce  to  those  of  [13]  if  the  F/s  are  chosen  in  the  form  of  rational 
polynomials  in  the  z-transform  variable  that  have  a  zero  at  z  =  0. 
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The  object  of  the  identification  process  is  to  adjust  the  para¬ 
s'  •  .  2 

meters  of  the  model  to  minimise  the  mean-squared  error,  E  -  ;|e:;  . 

An  estimate  of  the  process  transfer  function  is  found  fay  using  those 
minimizing  parameters  along  with  filter  transfer  functions,  F^,  in 
(3) .  Routinely,  the  parameters  that  minimize  the  error  norm  in  the 
approximate  model  are  used  in  (3)  even  though  it  is  the  transfer  func¬ 
tion  of  a  different  model.  The  result  is  an  estimate  of  the  process 
transfer  function  whose  error  depends  upon  tie  strength  of  the  noise 
process  "n". 

Why  bother  with  the  approximate  model  if  better  estimates  can  be 
had  using  the  actual  model?  This  question  is  quickly  answered  when 
one  attempts  to  estimate  the  parameters  of  the  actual  model.  Because 
the  model  output  and  the  parameters  are  interdependent,  the  parameter 
estimation  problem  is  highly  nonlinear,  and  must  be  solved  by  itera¬ 
tive  methods.  In  the  approximate  model,  the  parameter-invariant  pro¬ 
cess  output  replaces  the  model  output,  thereby  making  the  problem 
linear.  Hence,  the  approximate  model  is  widely  preferred. 

2.3  Estimation  of  the  Model  Parameters 

An  estimator  is  a  rule  for  choosing  the  model  parameter  values  in 
a  way  that  tends  to  minimize  the  error  between  the  process  output  and 
the  model  output.  Two  types  of  estimators  are  distinguished:  discrete- 
time  and  continuous-time. 

The  discrete  estimator  operates  on  sampled  data  and  produces  an 

estimate  of  the  pulse  transfer  function  of  the  process.  The  discrete 

T 

mean-square  error  is  defined  as  E  =  where  £  =  [e(0)  ...  e(M-D] 
is  a  column  vector  consisting  of  the  sampled  sequence. 

Likewise,  tne  model  output  sequence  and  filter  output  sequences 
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are  denoted  as  M- length  vectors  y^  and  y^i  for  i=l,...,n,  respec¬ 
tively.  The  output  of  the  approximate  model  can  be  expressed  as 


L n 


n 

Z 

i=l 


~ct.  y.  +g.u. 
1  1—1 


where 


n  = 


y, (0)  . . .  y  (0) 
n 

•  • 

y, (M-l)  ...  y'(M-l) 
i  n 


un (0)  ...  u  (0) 

■i  n 

*  • 

•  • 

u (M-l). ..  u  (M-l) 
1  n 


and  _8  is  the  vector  of  model  parameters 
8  =  [-aj...  -an81 

The  process  output  is  denoted  by  an  M- length  vector  The  error  can 
then  be  written  as 


n  e 


o 


It  can  be  shown  that  the  partial  derivatives  of  E  with  respect  to  the 
real  parts  of  the  parameters  are 


and 


8  E 

3  Rea . 
2 

3  E 

3  Reg. 
2 


2  Re [-e^y  ]  , 


a 

o 

2  Rete'u^J , 


a 


o 

for  j=l,...,n.  A  necessary  condition  that  a  set  of  parameters  must 
satisfy  in  order  to  minimize  E  is  that  the  partial  derivatives  listed 
above  vanish.  This  leads  to  the  set  of  normal  equations 

e'Vj  =  0, 

e'u.  =  0, 

-  “J 
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for  j=l,...,  n  which  can  be  written  in  matrix  form  as 


o'e  =  0 


or 


ne. 

=  0 


Solving  this  equation  for  _6  produces  the  discrete  least-squares 
estimator: 

=  [STa]-1  H'-a  v  . 

The  columns  of  a  consist  of  the  sampled  filter  outputs  of  the  approxi¬ 
mate  model  indicated  in  Figure  2.  The  filter  output  sequences  can  be 
computed  with  a  set  of  difference  equations  that  characterize  the  dis¬ 
crete  filters,  F^. 


If  H=2n,  (4)  reduces  to 


= 


e  »  a 


o  y 
o 


(5) 


If  M<2n,  ft'ft  is  singular  because  insufficient  data  are 
available  to  estimate  2n  parameters. 

The  continuous  estimator  operates  on  continuous  data.  The 
filters  within  the  model  are  continuous  in  this  case.  The  continuous 
mean-square  error  is  defined  as 


E  -  <e 


e*(t)  e(t)dt 


where  t^  and  t^  define  the  interval  on  which  the  waveform  exists.  The 
continuous  least-sauare  estimator  has  the  form: 


9  -  G"1  a  r 
—  2n  o  — 


(6a) 


where 


•••  ~ 3,  ■  •  •  3  ]  , 


(6b) 
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L=l<Y,y^  *••<?’>,)  <y’u]>  •••<y’un)]T*  (6c) 

<yi'7;>  -<V^  <Vyi>  •••<“„ 

<yryn>  •  •  • 

<fyl'ul>  •  •  • 

<yi;“r,>  •  •  •  <Un 

The  inner  product  is  defined  as 

/*t2 

=  J  ^  (t)  X2(t)  dt  .  (6e) 

C1 

The  continuous  estimator  results  from  the  minimization  of  the 
error  in  much  the  same  way  the  discrete  estimator  does,  and  therefore, 
the  deriviation  for  the  continuous  estimator  will  not  be  presented. 

The  parameters  of  the  actual  model  can  be  estimated  by  an  iter¬ 
ative  procedure  using  a  modification  of  the  estimators  described  above. 
This  procedure  involves  replacing  the  process  output  with  the  model 
output  as  the  feedback  source  within  the  model.  Specifically,  if 
y\^  replaces  y^  replaces  and  8^  replaces  in  (4)  and  (6), 
for  i  =  1,...,  n,  a  technique  of  iterative  improvement  results.  The 
parameters  of  the  iteration,  cu  and  are  estimated  in  terms  of 

the  most  recent  model  output,  y^  An  initial  estimate  of  the  orocess 

in 

is  required  to  start  this  procedure.  The  approximate  model  can  be  used 
to  provide  this  estimate.  Nothing  is  known  about  the  convergence 
characteristics  of  this  procedure.  However,  Steiglitz  and  Mc3ride  [13] 


do  mention  that,  in  general,  this  procedure  fails  to  converge  if  the 
initial  parameters  are  far  from  the  optimum  values. 


CHAPTER  III 


EXISTING  METHODS 


3.1  Introduction 

By  necessity  the  models  described  in  the  previous  chapter  are 
quite  general  since  they  must  serve  as  a  common  ground  for  the  two 
existing  methods  that  are  examined  in  this  chapter.  Prony's  method 
and  the  pencil-of-functions  method  are  shown  to  be  special  cases  of 
the  estimation  method  for  the  approximate  model.  Each  method  corres¬ 
ponds  to  a  particular  choice  of  filters  within  the  model. 

3.2  Prony's  Method 

In  1795  a  method  for  exponential  analysis  was  invented  by 
R.  Prony  [24].  The  method  was  implemented  by  hand  calculations. 

Prony's  method  has  since  found  numerous  applications  and  has  sometimes 
been  published  without  reference  to  its  inventor.  Van  Blaricum  [25] 
has  compiled  a  large  though  abridged  bibliography  of  Prony's  method. 

Van  Blaricum  and  Mittra  [8,9]  conducted  a  well  documented  investiga¬ 
tion  of  Prony's  method  and  proposed  some  useful  extensions  of  the 
technique.  Lager,  et  al.  [26]  have  suggested  a  "sliding-window  Prony" 
procedure  as  a  means  of  reducing  the  noise  sensitivity  of  the  method. 

Although  Prony's  method  is  computationally  efficient,  Dudley  [11] 
demonstrated  that  the  least-squares  version  [9]  can  produce  biased  pole 
values  that  differ  significantly  from  the  "best,,p°le  values  unless  the 
noise  component  of  the  waveform  is  small.  The  "best"  pole  values  are 
those  that  best  fit  the  waveform.  Prony's  method  results  from  the 
approximate  model.  It  follows  that  the  bias  is  nothing  more  than  the 
slight  error  that  is  incurred  by  the  use  of  the  approximate  model. 
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Prony's  method  results  from  the  approximate  model  under  the 

following  assumptions: 

1.  a  =1. 

n 

2-  F.  *  F.(z)  »  2 \ 

3.  The  model  input  is  an  impulse  at  k  =  0. 

4.  M  •  2n. 


The  assumption  of  «  *i  implies  that  the  parameters  are  normalized  rel- 
n 

ative  to  which  is  a  distinguishing  characteristic  of  Prony's  method 
The  other  methods  examined  in  this  work  are  normalized  relative  to 

aQ»  that  is*  aQ=i  is  assumed.  Since  M  =*  2n,  (5)  may  be  used  to 
estimate  the  parameters.  In  this  case. 


1  =  [-c^ - an_x  -1  Sx  3q]T  , 


y  5  [y(0)  •••  y(2n-l)r. 


and  ft  = 


y(l)  •••  y(n) 


v(n)  ♦ • «  y(2n-l) 
y(n+l)  y(2n) 


y(2n)  ••*  y(3n-l) 


I 


The  rank  of  si  is,  at  most,  n.  Hence,  with  some  rearranging  the 
following  equation  results: 


System  (7)  is  equivalent  to  Prony's  method  for  estimating  the  poles  of 

a  waveform.  The  poles  of  the  process  may  be  determined  from  the 

The  3^  must  be  determined  by  other  means.  See,  for  instance,  reference 

fllj.  A  straightforward  way  of  determining  the  A^  of  (1)  is  to  use 

the  poles  to  form  exponential  basis  functions  which  are  used  in  linear 

combination  to  form  a  least-squares  fit  to  the  data.  The  A^  are  the 

coefficients  of  this  linear  combination.  The  A may  also  be  found 

from  the  cr  and  8^,  once  these  parameters  are  known,  as  simply  the 

residues  of  H  . 

a 

Examples  of  the  performance  of  Prony's  method  under  a  variety  of 
conditions  are  found  in  sufficient  detail  elsewhere  [8,  11]  and,  for 
this  reason,  will  not  be  given  here. 

3.3  The  Peocil-of~Functions  Method 

The  pencil-of- functions  method  was  originally  proposed  by  Jain 
and  Gupta  [16]  before  Prony's  Hethod  became  well  known.  In  this 
section  the  pencil-of-functions  method  is  derived  as  a  special  case 
of  the  approximate  model  presented  in  Chapter  II  in  contrast  to  its 
usual  derivation  from  the  so-called  pencil-of-functions  concept  [19]. 
Additional  material  on  the  pencil-of-functions  method  is  contained 
in  references  [14,  15,  17],  A  discrete  version  of  the  method  is 
described  in  reference  [13], 

The  pencil-of-functions  method  results  from  the  approximate  model 
under  the  following  assumptions: 

1.  a  -  1. 

o 

2.  F.  =  F.  (s)  =  (1/s)1. 

i  i 

3.  The  matrix  G7n+1 ,  defined  in  what  follows,  is  singular. 

Note  chat  the  approximate  model  reduces  to  the  generalized  model  [1,3, 
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4]  under  the  first  assumption;  thus,  it  follows  that  the  pencil  of- 
functions  method  may  be  derived  from  the  generalized  model.  The  esti¬ 
mator  embodied  in  equation  (6)  may  be  used  to  estimate  the  transfer 
function  of  the  process,  and  the  poles  of  the  process  may  be  esti¬ 
mated  as  the  zeros  of  the  denominator  of  (3),  viz.. 


1  +  31  ur  **•  +  an 


r- 


or  Is11  +  a.  s11*1  +*-*+3=0.  (8) 

1  n 

The  may  be  found  using  the  estimator  of  (6).  3v  applying  Cramer’s 
rule  to  system  (6) ,  it  can  be  shown  that 

*  40iM00  (9: 

where  denotes  the  cofactor  formed  by  deleting  rcw  pt-1  and  column 
q+i  of  the  matrix. 


<™>  <y.y^>  ...  (y,y^  <r,»^  ...  (y,^ 

<h’P 


G2«l-  <7n'?> 


-(10) 


<V?> 


f’"3 


In  this  case  G„_ has  only  pure  real  elements.  If  G„  ,  is  also 
zirr±  Zn-H. 

singular,  the  following  relation  holds  among  the  cofactors: 


a  /  A 

ca  _  ^  l  ao 


This  relation  is  proved  in  the  Appendix.  The  ?encil-of-f unctions  method 
assumes  that  G9n+1  is  singular  and  makes  use  of  (8),  (9),  and  (II)  to 
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form  Che  following  characteristic  equation  for  the  process: 


It  may  be  shown  that,  in  general,  is  nonsingular  because  y  cannot 

be  expressed  as  a  linear  combination  of  the  y^  and  u^  under  the  fr  .low¬ 
ing  conditions  of  imperfect  modeling: 

1.  y  has  a  random  noise  component. 

2.  The  model  order,  n,  is  less  than  the  order  of 
the  process  or  waveform  being  modeled. 

3.  Both  1  and  2. 

If  any  of  the  above  conditions  apply  the  set  of  functions, 

s  =  {y,  V".  V  V">u.}  ' 

is  linearly  independent,  and  due  to  this  fact,  the  Gram  matrix  [27], 
^2n+l’  formec*  from  those  functions  of  S  is  nonsingular.  Otherwise, 
for  conditions  of  perfect  modeling,  the  functions  of  S  are  linearly 
dependent,  G„  ,,  is  singular,  and  (11)  holds. 

Thus,  if  the  process  can  be  modeled  perfectly,  that  is,  with  no 
error,  then  (12)  is  equivalent  to  (8).  However,  if  the  process  cannot 
be  modeled  perfectly,  (11)  no  longer  holds  and  (12)  becomes  less 
accurate  than  (8)  for  estimating  the  parameters  of  the  approximate 
model.  Since  the  para^-'cer  values  of  the  approximate  model  are  trans¬ 
planted  into  the  true  model,  it  is  not  clear  if  the  a.  or  the  /S77 75 

x  xx  oo 

more  accurately  portray  the  process  after  this  transplant.  In  the 

next  section  numerical  evidence  is  presented  that  indicates  that  the 

provide  better  estimates.  However,  since  the  a.  and  the  /  A . .  /  A  are 

1  11  oo 

computed  in  extremely  different  ways,  the  results  may  not  indicate  the 
true  situation  due  to  inaccuracies  in  computation. 
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Aside  from  the  problem  mentioned  above  the  pencil-of-functions 
method  has  two  other  difficulties  of  which  only  one  has  a  simple 
remedy. 

The  difficulty  that  has  a  simple  remedy  is  related  to  the  use  of 
integrators  for  the  filters  within  the  system  model.  The  continuous¬ 
time  integrators  cannot  be  implemented  exactly  by  any  algorithm.  When 
approximate  integrators  are  cascaded  as  they  are  in  the  pencil-of- 
functions  method,  large  errors  accumulate  quickly  and  the  intended 
result  destroyed  after  a  number  of  integrations.*  This  difficulty  can 
be  resolved  by  cascading  discrete  integrators  to  obtain  filters  with 
pulse  transfer  functions  given  by 


>.  ■  •  A>‘  •  M* 


which  can  be  implemented  on  a  digital  computer  with  no  error  by  using 
the  difference  equations: 

yi(k)  *  yi(k-l)  +  y^OO 


and  u^,(k)  =  u^Ck-l)  +  u.^Ck) 


The  variable  Z  is  defined  as 


and  z  is  the  z-transform  variable.  When  discrete  integrators  are  used 
the  discrete  pencil-of-functions  method  results.  The  poles  of  the  pulse 
transfer  function  of  the  process  may  be  estimated  as 


*This  feature  of  the  pencil-of-functions  method  has  been  consistently 
observed  in  numerical  implementations.  This  encroachment  of  systematic 
error  is  inconsistent  with  generally  accepted  interpretation  of  numeri¬ 
cal  integration  processes  and  has  not  been  explained  analytically,  to 
date. 
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where  Z.  is  the  ith  zero  of 

i 


;th 
,n-l 


iZn  +  a  ZQ  1  +  •  *  *  +  ot  =0  (13) 

1  n 

and  the  o  are  found  by  use  of  equation  (4) . 

If  the  continuous  inner  product  (6e)  is  replaced  with  the  discrete 
inner  product. 


<£l,x\?  =  X*  (k)  X2'  (k)  .  (14) 

*  k=0 

It  can  be  shown  that  rD  »  Q'z  and  G2n  3  where  and  r.  are  formed 

with  the  discrete  inner  product  (14)  in  the  same  way  that  G2n  and  r  are 
formed  with  the  continuous  inner  product  (6e) .  Then, 


a 


i 


(15) 


where  A  denotes  the  cofactor  formed  by  deleting  row  p+1  and  column 
q+1  of  G°  . , ,  and  G?  ..  is  the  Gram  matrix  formed  with  the  discrete 
inner  product  in  the  same  way  G2n+1  ^ormed  with  the  continuous 
inner  product.  Equation  (15)  follows  from  Cramer's  rule  in  the  same 
way  equation  (9)  does.  As  before  the  relation, 

TT 


ad 

SI 

A° 

PP 


(16) 


holds  if  G, 


!n+l 


is  singular.  Combining  (13),  (15),  and  (16)  yields  the 


characteristic  equation  of  the  system: 
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Equation  (17)  is  the  equivalent  of  (12)  for  the  discrete  method.  The 
choice  of  positive  radicals  in  (17)  follows  from  arguments  counterpart 

to  those  put  forth  by  Jain  [19]  for  the  continuous  case. 

Estimates  of  the  poles  of  the  Laplace  transfer  function  of  the  pro¬ 
cess,  s  ^ ,  are  related  to  the  zeros  of  (13)  and  (17)  by 
*n(l-Z.) 

Si=  T  ' 

In  the  next  section  the  discrete  version  of  the  pencil-of-functions 
method  is  tested  on  noisy  data  and  the  performance  of  equation  (13)  is 
compared  to  that  of  (17)  in  order  to  demonstrate  the  superiority  of  (13) 
in  estimating  the  system  poles. 

Now  attention  is  turned  to  the  second  difficulty  with  the  method 
that  does  not  have  a  simple  remedy.  This  difficulty  is  related  to 
the  attenuation  of  the  higher  frequency  modes  of  the  process  output  by 
the  repeated  integrations  applied  to  the  output  waveform.  It  can  be 
verified  that  an  integrator  is  simply  a  first  order  filter  whose 
Laplace  transfer  function  has  a  pole  at  the  origin  in  the  s-plane. 

Such  a  filter  tends  to  suppress  the  higher  frequencies  present  at  its 
input.  The  higher  frequency  suppression  phenomenon  is  illustrated  in 
Figure  3.  Normally,  when  an  exponential  function  is  integrated  re¬ 
peatedly,  components  consisting  of  powers  of  time  exist  in  the  higher 
integrals  as  well  as  the  original  exponential  function  components. 

In  Figure  3  the  components  of  powers  of  time  have  been  subtracted  from 
the  integrated  waveforms  in  order  to  make  the  attenuation  of  the  higher 
modes  more  evident.  The  first  waveform  is  a  hypothetical  waveform  pro¬ 
vided  for  analysis.  The  waveforms  that  follow  are  the  integrals  of 
increasing  order  of  the  first  waveform  and  display  the  increasing 
dominance  of  the  fundamental  mode  or  the  mode  of  lowest  frequency. 
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Further  integrations  yield  nearly  identical  waveforms.  The  integrated 


waveforms  tend  to  become  linearly  dependent  at  higher  model  orders. 


D 


The  Gram  matrix,  G2q  or  G^,  then  tends  to  singularity  and  the  method 


becomes  unstable.  The  suppression  phenomenon  occurs  in  both  the 
discrete  and  continuous  methods.  In  fact,  even  for  very  modest  model 
orders,  the  method  can  become  numerically  ill-conditioned  to  a  degree 


that  special  care  must  be  taken  to  assure  accurate  inversion  of  G^^  or 


G„  .  The  examples  of  the  next  section  illustrate  the  increasing  ill- 
zn 


condition  with  increasing  model  order  by  computing  the  condition  number 


o£  G2n • 


3.4  Numerical  Examples  of  the  Pencil-of-Functions  Method 

In  the  first  example,  impulse  input  to  the  system  is  assumed  and 
the  discrete  method  is  applied  in  the  analysis  of  a  waveform  consisting 
of  50  samples  defined  by: 


y(k)  =  ^  A^eSjkT  +  n(k) 


k  »  0,1, •••,49 


where  =  1, 


sx  -  -1  +  jlO, 


a2  =  1, 


s7  =  -1  -  jlO, 


A^  3  1, 


=  -1.5  +  j30. 


\  ~  1» 


s4  =  -1.5  -  j30. 


n(k)  is  a  Gaussian  distributed  white  noise  sequence  of  50  samples,  and 
T  =  1/49.  The  signal-to-noise  ratio  is  30  dB,  where  signal-to-noise 
ratio  (SNR)  is  defined  by: 


SNR  (dB)  =  20  log1Q  ~ 


(18) 


i  m 


5  , 


max 

and  a=i. 


(maximum  magnitude  of 
uncorrupted  waveform) 


A.esjkT 

3-1  : 

o  =  standard  deviation  of  noise. 

Four  poles  are  requested.  Thirty  Monte  Carlo  runs  are  made  where  each 
run  uses  a  different  noise  sequence.  Figure  4  shows  the  resulting 

s-plane  pole  estima. for  all  Monte  Carlo  runs  overlaid  onto  a  common 

plot.  In  Figure  4(a),  the  poles  are  estimated  with  the  roots  of  equa¬ 
tion  (17),  and  in  Figure  4(b),  the  poles  are  estimated  with  the  roots 
of  equation  (13) .  The  reader  is  reminded  that  the  residues  of  both 
poles  are  of  unit  amplitude  thus  leading  to  the  conclusion  that  the 
apparent  deterioration  of  pole  accuracy  with  increasing  frequency  is  an 
artifice  of  the  process. 

In  the  second  example  the  same  analysis  is  carried  out  for  a 
higher  model  order.  The  waveform  is  given  by: 

6 

v(k)  =  A.eSJkT  +  n(k)  ,  k  =  0,1, •••,49  , 

where  A^  =1,  s,  =  -2  +  j50, 

A6  =1,  s&  =  -2  -  j5Q, 

and  all  other  details  remain  unchanged  from  the  first  example.  Six  poles 
are  requested.  Figure  5  shows  the  resulting  S-plane  pole  estimates.  In 
Figure  5(a)  the  poles  are  estimated  with  the  roots  of  (17),  and  in  Fig¬ 
ure  5(b)  the  poles  are  estimated  with  the  roots  of  (13). 

The  condition  number  of  is  averaged  over  the  thirty  Monte  Carlo 
runs  and  displayed  in  Figures  4  and  5.  The  condition  number  is  defined 


(a)  estimates  using  roots  of  (17) 

(b)  estimates  using  roots  of  (13) 


(a)  estimates  using  roots  of  (17), 

(b)  estimates  using  roots  of  (13). 
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Condition  number  =  |  [G,?  j  j  j  |  [G.?  ]  |  | 


where  | j  •  | |  denotes  the  Chebyshev  matrix  norm  defined  by: 

N 

i  iai  i  V  |a  i 

'  l<i<N  ij  1 

til  til 

and  a.,  is  the  element  of  the  i  row  and  j  column  of  the  N-dimensional 
ij 

matrix  A.  The  larger  the  condition  number,  the  more  ill-conditioned  the 
matrix  A. 

The  results  of  these  examples  demonstrate  that  better  pole  esti¬ 
mates  are  obtained  by  using  (13)  instead  of  (17) .  The  increasing  con¬ 
dition  number  from  the  first  to  the  second  example  indicate  that  the 
method  is  becoming  rapidly  ill-conditioned  as  the  order  of  the  method  is 
increased.  The  results  of  a  third  example  for  which  an  eighth  order 
method  was  used  to  analyze  a  noise  corrupted  waveform  composed  of  eight 
exponentia]  components  are  not  presented  because  the  method  became  so 
ill-conditioned  that  the  intended  result  was  completely  destroyed. 

It  was  found  that  the  pencil-of-functions  method  becomes  ill- 
conditioned  for  the  analysis  of  waveforms  containing  many  poles.  Even 
though  the  results  of  these  examples  do  not  apply  to  the  continuous 
version  of  the  method  directly,  it  is  known  that  in  the  limit  as  T 
approaches  zero  the  discrete  method  approaches  the  continuous  method, 
and  hence,  there  is  a  measure  of  similarity  in  the  two  methods.  This 
measure  of  similarity  is  felt  to  be  sufficient  to  claim  that  the  con¬ 
tinuous  version  of  the  method  becomes  ill-conditioned  at  higher  order. 


Chapter  IV 


THE  ADAPTIVE  METHOD 


4.1  Introduction 

The  adaptive  method  is  an  iterative  method  which  is  quite  tolerant 
of  noise  in  the  waveform  undergoing  exponential  analysis.  The  original 
idea  for  this  method  was  inspired  by  the  pencil-of-functions  concept 
due  to  Jain  [17].  In  this  chapter,  the  adaptive  method  is  derived,  not 
from  the  pencil-of-functions  concept  as  it  was  originally,  but  rather 
from  the  parametric  models  cf  Chapter  II.  Formulating  the  method  in 
this  way  is  not  only  simpler  but  yields  greater  insight  into  the 
method's  nature. 

4.2  An  Adaptive  Filtering  Scheme 

The  adaptive  method  results  from  the  identification  scheme  of 
Chapter  II  under  the  following  assumptions: 


^i  z-z. 

i 

3.  The  model  input  is  a  unit  sample  at  k  *  0  (discrete 
impulse) . 

The  unique  feature  of  the  method  is  that  the  filter  poles,  z^,  may  be 
adjusted  to  any  value  in  the  z-plane.  An  adaptive  technique  for  ad¬ 
justing  the  filters  consists  of  first  initializing  the  filters  to 
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arbitrary  values  in  the  z-plane  and  repeating  the  following  steps: 


1.  Find  an  estimate  of  the  process  transfer  function 
using  the  current  filters  in  the  model. 

2.  Set  each  filter  pole  to  one  pole  of  the  estimated 
transfer  function. 

This  process  is  repeated  until  the  approach  zero. 

It  is  shown  in  the  next  section  that  the  poles  of  the  process 
can  be  estimated  during  the  course  of  iteration  by 
z. 


2  .  = 


'i  l  +  o. 

i 


(19) 


removing  the  need  for  finding  the  roots  of  a  polynomial.  The  filter 
poles  are  updated  to  z.  on  each  iteration.  When  the  o.  approach  zero, 
the  pole  updating  ceases  and  the  method  converges.  At  each  iteration, 
the  a^  are  found  by  using  equation  (4) .  At  convergence  the  3-plane 
poles,  s^,  can  be  obtained  from  the  filter  poles  by 
2n  z. 


s  = 


(TO) 


and  the  A.  =  B.. 

i  i 


4. 3  Estimation  of  Process  Poles  from  Transfer  Function  Parameters 
As  usual,  the  poles  of  the  process  transfer  function  can  be 
estimated  as  the  poles  of  (3).  However,  to  estimate  the  poles,  a 
polynomial  in  z  is  needed.  To  find  this  polynomial  the  numerator  and 
denominator  of  (3)  are  multiplied  by  the  product  of  all  filter  denomi¬ 
nators,  that  is,  by  ^  .  The  required  polynomial  is  the  de- 

i=l  1 

nominator  of  the  expression  that  results  after  the  multiplications 


-31- 


indicated  above  are  carried  out  and  terns  of  equal  degree  have  been 
combined  in  the  denominator.  Although  it  can  be  done,  the  formulation 
of  the  polynomial  coefficients  on  a  digital  computer  requires  some 
small  computational  burden.  But  another  method  exists  for  estimating 
the  poles  that  bypasses  this  complication  as  well  as  the  necessity  for 
finding  the  roots  of  a  polynomial.  For  the  sake  of  discussion,  assume 
that  the  filter  poles  have  been  adjusted  to  nearly  coincide  with  the 
process  poles  and  examine  the  behavior  of  the  (i+1) term  of  the  de- 
nominator  of  (3)  when  the  variable  2  approaches  the  value  of  the  i 
process  pole.  It  is  observed  that  the  (i+l)C^  term  can  become  arbitrar¬ 
ily  large  if  the  ittl  filter  pole  is  arbitrarily  close  to  the  iCil  process 
pole.  It  is  then  possible  to  write  a  simplified  expression  for  the 

denominator  of  (3)  which  is  approximately  equivalent  when  2  is  equal 
cil 

to  the  1  process  pole,  viz. , 


1 


V. 

2-2. 


(21) 


where  all  other  terms  other  than  the  (i+l)th  term  in  the  denominator 
of  (3)  are  negligibly  small.  If 
a.  z. 

1  +  —  ■  ■  =  0  (22) 

£r2i 

where  2^  is  the  value  of  z  for  which  the  equality  holds,  then  z^  is  a 
reasonable  estimate  of  the  iC^  process  pole.  Equation  (19)  is  obtained 
by  solving  (22)  for  z^.  The  approximations  of  (21)  and  (22)  follow 
from  observing  that  near  convergence  of  the  z^  all  of  the  except  aQ, 

which  is  fixed  as  unity,  vanish.  Therefore,  the  aQ  term  dominates 
along  with  the  (i-HL)  term  as  argued  above. 
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4.4  Xumericai  Examples  of  the  Adaptive  Method 

The  first  example  of  the  adaptive  method  uses  a  simple  first-order 

model  to  analyze  a  sequence  consisting  of  three  samples.  Because 

of  the  example’s  simplicity,  hand  calculations  and  concrete  numbers 

c3n  be  displayed. 

The  sequence  to  be  analyzed  is 

{y(k),  k-0,1,2}  =  {2,1,2}  . 

It  can  be  shown  that 

y  (k)  *  4  (a  constant  sequence) 
m  J 

minimizes  the  mean-squared  error  between  y  and  y  if  the  problem  is 

m 

restricted  to  a  first-order  solution.  Therefore,  the  expected  results 
of  this  analysis  after  the  adaptive  method  converges  are:  ~  and 

z^  =  1.  Let  the  initial  filter  pole,  z^,  be  set  to  one.  One  might 
expect  that  the  method  would  be  convergent  immediately  in  this  case. 
However,  the  method,  in  fact,  will  not  converge  at  this  value  of  z1 
due  to  the  error  introduced  by  the  use  of  the  approximate  model. 

The  convergent  value  of  z^  should  be  slightly  perturbed  from  the 
expected  value. 

The  estimator  of  equation  (4)  is  applied  to  obtain  estimates  of 
and  8^-  In  this  case. 
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2  1 
3  1 
5  1 


and  the  system  that  must  be  solved  is 


The  estimates  are  a^=-l/14=-. 07143  and  8^=20 /14~1. 429.  The  estimate  of 

the  process  pole  is 

Z1  1  14  -  .  n,_ 

Z1  l+aL  "  1-1/14  13  '  1,077 ‘ 


The  estimates  above  constitute  the  results  of  the  first  iteration. 
For  the  second  iteration  the  filter  pole  is  updated  to  1.077  and  the 
same  procedure  is  repeated  for  iterations  two  through  five.  The 
results  for  all  iterations  are  displayed  in  Table  1. 

Next,  the  iterative  scheme  for  the  actual  model  set  forth  in 
Chapter  II  is  simultaneously  combined  with  the  adaptive  filtering 
scheme  and  is  applied  to  correct  the  biased  value  just  obtained  with 
the  approximate  model.  An  initial  estimate  of  the  process  transfer 
function  is  required  to  start  this  estimation  procedure.  The  para¬ 
meters  obtained  with  the  approximate  model  at  convergence  are  used 

to  form  this  initial  estimate.  The  initial  pole  is  1.077.  The 

th 

most  recent  model  output  sequence  is  computed  for  the  n  order  ^ase 

by  using  the  following  difference  equations  in  the  order  indicated: 

Ui(k)  =  u(k)  +  2.u.(k-l), 
i  ii 


v 

'  m 


(k) 


J(8iu.(k)-Wmi(k-1>) 


n 

1  +  :  a. 
i-1 1 
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Table  1  Estimates  using  the  approximate  model. 


Iteration 

number 

-°1 

S1 

21 

1 

7.143  x 

l(f2 

1.429 

1.077 

2 

2.711  x 

io"4 

1.539 

1.077 

3 

2.279  x 

l(f6 

1.539 

1.077 

4 

-4.245  x 

io-7 

1.539 

1.077 

5 

2.954  x 

io" 7 

1.539 

1.077 

Table  2  Sequences  on  the  first  iteration  of  the  estimation  scheme 
for  the  actual  model. 


k 

ux(k) 

\<k> 

0 

1 

1.539 

1.539 

1 

1.077 

1.658 

3.316 

2 

1.160 

1.785 

5.  356 

! 


y_,  00  =  y  (k)  +  z  y  . (k-1). 
mi  m  x  mi 

cti 

These  equations  are  valid  for  the  n  order  model.  For  the  first 
order  (n*l,  i=l)  model  that  is  considered  here,  the  equations  are 
u^Ck)  =  u(k)  +  z^u^Ck-l), 


y  (k)  = 

m 


1  +  a , 


yml(k)=  ym(10+  Vml0^’ 

and  are  used  to  compute  the  results  shown  in  Table  2.  From  the  results 
of  Table  2, 


Q 


1.539  1.000 

3.316  1.077 
5.356  1.160 


and  the  system  which  must  be  solved  is 


42.051  11.323  -a  _  17. 

11.323  3.506  2^  5. 


106 

397 


The  estimates  are  =  +5.972x10“^  and  3^  5  1.732.  The  estimate  of 
the  process  pole  is 

- _ ZJ 

Zl  ~ 


1.077 


1  +  a. 


1  +  .05973 


=  1.017. 


The  estimates  just  found  constitute  the  results  of  the  first  itera¬ 
tion  using  the  estimation  scheme  for  the  actual  model.  The  filter 
pole  is  updated  to  1.017  and  the  same  procedure  is  repeated  several 
more  times  to  obtain  the  results  shown  in  Taole  3. 


The  results  in  Table  3  indicate  that  the  estimation  procedure  for 
the  actual  model  has  indeed  converged  to  the  expected  parameter  values 


i 
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Table  3  Estimates  u~  •  he  actual  model. 


Iteration 

number 

“al 

B1 

21 

1 

-5.973  x 

10~2 

1.732 

1.017 

2 

-2.053  x 

10~2 

1.705 

.996 

3 

1.163  x 

10-3 

1.669 

.997 

4 

2.349  x 

10~3 

1.663 

1.000 

5 

5.878  x 

io"4 

1.665 

1.000 

6 

-5.633  x 

10~5 

1.667 

1.000 

7 

-7.576  x 

io'5 

1.667 

1.000 

8 

-1.697  x 

io"5 

1.667 

1.000 

9 

2.032  x 

10~6 

1.667 

1.000 

10 

2.346  x 

io'6 

1.667 

1.000 

11 

9.632  x 

io'7 

1.667 

1.000 

12 

-2.289  x 

io~7 

1.667 

1.000 

13 

-1.717  x 

10~7 

1.667 

1.000 

14 

2.861  x 

io-8 

1.667 

1.000 

i  * 
i  : 
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whereas  the  results  in  Table  1  clearly  show  that  the  procedure  for  the 
approximate  model  does  not.  However,  the  magnitude  of  the  error  that 
the  use  of  the  approximate  model  introduces  in  the  estimated  pole 
value  is  relatively  small. 

The  results  of  the  next  example  provide  further  evidence  that, 
indeed,  the  error  in  the  parameter  estimates  are  small  when  adaptive 
filters  are  used  in  the  approximate  model. 

To  illustrate  the  performance  of  the  adaptive  method  on  measured 
data,  the  method  was  applied  in  the  extraction  of  the  natural  resonance 
of  the  electrical  transient  response  of  a  thin  cylinder  in  free  space. 
Figure  6  shows  the  responses  of  the  cylinder  at  five  points  along  its 
length  which  were  measured  by  techniques  described  in  reference  [29]. 

The  cylinder  was  excited  by  a  500  picosecond  burst  of  radiation  that  is 
normally  incident  to  its  axis.  The  cylinder  of  60  centimeters  long  and 
approximately  1  centimeter  in  diameter.  The  waveform  consists  of  512 
samples  and  has  a  time  step  of  .9775x10  11  seconds.  The  first  109  sam¬ 
ples  are  ignored  since  the  forcing  wavefront  impinging  on  the  cylinder 
corrupts  this  portion  of  the  natural  response.  The  preprocessing  incre¬ 
ment  is  10  samples.  This  means  that  every  10  adjacent  samples  are 
averaged  to  form  one  sample  of  the  preprocessed  waveform  beginning  at 
sample  110.  The  preprocessed  waveform  then  has  the  integer  portion  of 
(512-109) /10  or  40  samples  with  a  time  step  of  (10)  •  (.0775x10  ^)  = 
.9775xl0“10  seconds.  The  estimation  procedures  for  both  the  approx¬ 
imate  model  and  the  actual  model  were  applied  to  the  preprocessed 
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waveforms  and  the  results  for  each  case  is  shown  in  Figure  7.  The 
extracted  poles  from  all  responses  are  shown  overlaid  in  Figure  7. 
Although  the  data  appear  to  be  free  of  noise,  the  signal  to  noise 
ratio  is  estimated  to  be  between  15  dB  and  20  dB  by  observing  the 
model  error.  The  noise  is  thought  to  arise  from  such  phenomenon 
as  reflections  on  the  transient  range  and  in  the  probe  cabling  al¬ 
though  no  tests  were  made  to  confirm  this.  It  should  be  pointed 
out  that  the  poles  obtained  with  the  actual  model  actually  minimize 
the  mean-squared  error,  and  therefore  provide  the  best  least-squares 
fit  to  the  preprocessed  responses.  If  better  parameter  estimates  are 
to  be  found,  more  must  be  known  about  the  noise  that  corrupts  the 
waveforms.  It  is  interesting  to  note  that  the  results  obtained  with 
the  approximate  model  and  those  obtained  with  the  actual  model  are 
almost  identical.  One  has  to  strain  to  see  the  difference.  Many 
other  cases  have  been  studied  whose  results  have  confirmed  the  trend 
of  nearly  identical  pole  estimates.  In  many  practical  cases,  one  may 
choose  to  use  the  results  from  the  approximate  model  without  bothering 
to  refine  those  estimates  further  by  using  the  estimation  procedure 
for  the  actual  model.  This  may  be  a  wise  choice,  particularly  in  view 
of  the  fact  that  the  adaptive  estimation  procedure  for  the  actual 
model  does  not  converge  in  many  cases  where  the  adaptive  procedure 
for  the  approximate  model  does  converge. 

In  the  next  several  examples,  the  analysis  of  numerically  generated 
transient  data  is  examined.  The  transient  data  were  obtained  using 
the  time  domain  computer  code  TWTD  [30] .  The  structure  modeled  by 
TWTD  was  a  thin  cylindrical  scatterer.  This  structure  was  chosen 


centimeter  cylinder 
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because  of  the  availability  of  comparative  results  obtained  by  Tesche 
[31] .  The  true  pole  values  for  the  first  six  odd  harmonic  modes 
of  the  cylinder  obtained  by  Tesche  with  an  integral  equation  technique 
are  listed  in  Table  4.  These  values  are  the  theoretical  values  which 
are  used  for  comparison  in  the  examples  that  follow. 

The  first  three  examples  with  the  TWTD  data  illustrate  the  results 
that  can  be  obtained  with  the  adaptive  method  at  three  different 
signal-to-noise  ratios  (SNR):  25  dB,  20  dB,  and  15  dB.  Figure  8 
illustrates  a  noise  contaminated  TWTD  waveform  consisting  of  256  samples 
with  a  15  dB  SNR  where  SNR  is  defined  in  the  spirit  of  equation  (8) . 

The  time  step  is  0.68020  x  10  ^  seconds.  The  waveform  is  that  of  the 
current  versus  time  history  at  the  center  of  a  one-meter  cylinder  which 
is  excited  by  a  voltage  source  offset  approximately  15  centimeters  from 
the  center.  The  voltage  source  has  a  Gaussiaa-pulse  time  history  where 
pulse  width  referred  to  the  1/e  level  of  the  function  is  approximately  .6 
nanoseconds.  Although  the  cylinder  has  an  infinite  number  of  modes, 
only  the  first  four  low  frequency  modes  dominate  the  response  since 
the  excitation  is  band  limited.  Figure  9  shows  the  results  of  the 
estimation  procedures  for  both  the  actual  model  and  the  approximate 
model  for  five  Monte  Carlo  runs  where  each  run  uses  a  different 
noise  sequence.  The  noise  is  Gaussian  distributed  and  uncorrelated. 

The  first  60  samples  are  ignored  to  exclude  the  influence  of  the  driv¬ 
ing  voltage.  The  waveforms  are  preprocessed  before  being  analyzed 
with  a  preprocessing  increment  of  5  samples  beginning  at  sample  70. 


Table  4 


Cylinder  pole  values  predicted  by  Tesche  [31] 


sL/c 


Real 

Imag 

-0.082 

0.926 

-0.147 

2.874 

-0.188 

4.835 

-0.220 

6.800 

-0.247 

8.767 

-0.270 

10.733 

Figures  10  and  11  show  the  pole  estimates  for  the  same  data  with  SNR's 
of  20  and  25  dB,  respectively.  Tables  5  and  6  tabulate  the  pole 
values  corresponding  to  Figure  9;  Tables  7  and  8  correspond  to  Figure 
10;  and  Tables  9  and  10  correspond  to  Figure  11. 

Several  things  should  be  pointed  out  about  these  results.  First, 
even  though  there  are  four  dominate  inodes  in  the  data  only  three  pole 
pairs  are  requested  for  the  15  and  20  dB  cases.  Only  three  pole  pairs 
are  requested  because  the  eighth  order  method  often  diverges.  In 
several  of  the  Monte  Carlo  runs  not  all  of  the  poles  are  in  conjugate 
pairs.  However,  most  of  the  unmatched  poles  had  values  close  to  those 
reported  by  Tesche.  The  runs  which  did  not  converge  are  so  labeled. 
Nonconvergent  runs  usually  produce  pole  estimates  that  differ  dras¬ 
tically  from  Tesche' s  results.  Hence,  a  reasonable  procedure  to 
handle  the  case  where  the  method  converges  but  yields  unmatched  poles 
might  be  to  simply  assume  that  each  unmatched  pole  possesses  a  con¬ 
jugate  companion  pole.  In  fact,  a  promising  procedure  to  eliminate 
this  problem  would  be  to  increase  the  model  order  once  for  each  un¬ 
matched  pole  at  convergence,  introduce  the  appropriate  conjugate  pole, 
and  continue  the  iteration  until,  hopefully,  convergence  is  achieved 
with  all  poles  occurring  in  conjugate  pairs. 

For  the  25  dB  case  four  pole  pairs  are  requested  and  four  con¬ 
jugate  pole  pairs  converged  in  each  case  except  one  which  diverged. 

The  results  for  the  offset-driven  data  indicate  that  the  adaptive 
method  is  able  to  provide  useful  results  even  in  noise  levels  of 
around  15  dB  and  even  if  the  convergent  poles  are  not  in  conjugate 
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In  che  next  two  examples  the  adaptive  method  is  applied  in  the 
analysis  of  che  noise  contaminated  TNTD  data,  the  current  at  the  center 
of  a  one-meter  cylinder.  The  exciting  voltage  generator  has  a  Gaussian 
time  history  with  a  pulse  width  at  the  1/e  level  of  approximately  .3 
nanoseconds.  Hence,  che  center-driven  waveform  contains  components 
which  are  higher  in  frequency  than  the  offset-driven  waveform.  Figure 
12  displays  the  noise-contaminated  center-driven  waveform  with  a  20  dB 
SNR.  The  noise  is  Gaussian  distributed  and  uncorrelated.  The  first  69 
samples  are  ignored  to  exclude  the  influence  of  the  driving  voltage. 

The  waveforms  are  preprocessed  before  being  analyzed  with  a  preprocess¬ 
ing  increment  of  5  samples,  beginning  at  sample  70.  Figures  13  and  14 
show  the  pole  estimates  for  the  center-driven  TWTD  data  with  SNR's  of 
20  and  25  dB,  respectively.  Tables  11  and  12  tabulate  the  pole  values 
corresponding  to  Figure  13;  and  Tables  13  and  14  correspond  to  Figure 
14. 

The  results  for  the  center-driven  data  display  the  sensitivity 
of  the  higher  frequency  modes  to  noise.  This  sensitivity  is  probably 
due  to  a  lower  relative  degree  of  coupling  to  the  higher  frequency 
modes  for  the  center-driven  case.  That  is,  the  higher  frequency  modes 
are  relatively  weak  for  the  center-driven  case  and  are  more  easily 
corrupted  by  the  noise.  For  the  data  with  a  25  dB  SNR  five  pole  pairs 
are  extracted  with  success.  As  usual,  a  few  Monte  Carlo  runs  diverged. 


Noise -contaminated  offset-driven  TWT1)  waveform,  SNU“lf>  .«H.  The  uncontaminated 
waveform  is  plotted  under  the  noisy  waveform  for  comparison. 
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Table  10  Poles  for  the  offset-driven  TWTD  data,  SNR  =  25  dB.  Resul 
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The  preprocessing  of  the  waveforms  was  done  primarily  to  reduce 
the  number  of  samples  that  the  method  is  required  to  process  in  order 
to  improve  efficiency.  The  preprocessing  also  reduces  the  noise  level 
of  the  waveform  although  the  information  content  of  the  waveforms  re¬ 
mains  the  same.  The  performance  of  any  method  should  be  judged  by  how 
well  it  uses  available  information  instead  of  by  how  large  a  SNR  it  can 
tolerate. 

Attempts  to  analyze  waveforms  consisting  of  highly  damped  expo¬ 
nential  components,  such  as  the  transient  responses  of  a  sphere,  were 
not  successful.  The  adaptive  method  does  not  converge  for  such  wave¬ 
forms  which  display  double  pole  characteristics,  that  is,  waveforms 
wx  '  components  of  the  form  t  exp(st).  Slight  modifications  to  the 
adaptive  method  might  allow  the  analysis  of  such  waveforms.  The  de¬ 
scription  of  these  modifications  will  have  to  wait  until  further  study 
is  completed  in  this  direction. 


Chapter  V 


CONCLUSIONS  AND  FUTURE  DIRECTIONS 


This  document  has  related  three  identification  methods,  namely 
the  pencil-of-f unction  method,  Prony's  method,  and  the  adaptive 
method  to  the  general  identification  scheme  presented  in  Chapter  II. 
Each  method  results,  fundamentally,  from  a  certain  choice  for  the  fil¬ 
ters  in  the  general  identification  model.  A  meaningful  way  to  sum¬ 
marize  the  relation  between  the  three  methods  is  to  plot  the  transfer 
functions  of  their  respective  filters  as  shown  in  Figure  15.  It  should 
be  pointed  out  that  the  filters  of  Prony's  method  and  the  pencil-to- 
functions  method  are  cascaded,  while  the  filters  of  the  adaptive  method 
are  not.  The  poles  of  the  discrete  filters  have  been  mapped  into  the 
S-plane  with  the  mapping  defined  by:  z  =  exp(sT).  The  filters  of 
Prony's  method  treat  all  frequencies  in  an  equal  manner;  the  filters 
of  the  pencil-of-functions  method  are  preferential  to  the  low  fre¬ 
quencies;  those  of  the  adaptive  method  possess  passbands  which 
are  adjustable. 

The  adaptive  method  is  a  new  method  which,  in  many  cases,  pro¬ 
vides  excellent  poles  estimates  under  difficult  conditions.  The 
method  is  unique  in  that  a  solution  to  a  polynomial  is  not  required 
to  find  estimates  of  the  process  poles.  The  method,  in  effect, 
"swallows"  the  polynomial  solver  in  its  own  iterative  pole-searching 
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scheme.  The  adaptive  method  is  thought  to  be  closely  related  to  the 
method  of  Steiglitz  and  McBride  [19]  since  both  models  filter  the  in¬ 
put  and  output  records  with  filters  whose  poles  are  the  most  recent 
estimates  of  the  process  poles.  However,  the  method  of  Steiglitz  and 
McBride  does  not  provide  pole  estimates  but  only  estimates  of  the 
process  transfer  function  during  the  course  of  iteration.  The  primary 
utility  of  the  adaptive  method,  in  the  author's  opinion,  lies  in  re¬ 
finement  of  predetermined  poles. 

The  pencil-of-function  method  was  found  to  be  ill-conditioned  for 

the  identification  of  high  order  processes.  Although  the  method  can 

provide  superior  estimates  of  low-order  processes.  It  is  shown  that 

use  of  the  >/a.  ./A  in  place  of  the  o.  was  less  accurate  in  general, 

11  oo  i 

than  using  equation  (4)  or  (6)  to  estimate  the  for  the  approximate 

model.  Since  the  parameter  values  themselves  are  transplanted  into 

the  actual  model,  it  is  not  clear  if  the  ./A  or  the  a.  more 

n  oo  1 

accurately  portray  the  process  after  this  transplant.  But  there  is  no 

reason  to  believe  the  /aTT/A  provide  a  better  estimate  after  the 

11  oo 

transplant  and  every  reason  to  believe  that  they  do  not. 

Dudley  [ll]  indicated  that  the  noise  sensitivity  of  least-squares 
Prony's  method  was  due  to  parameter  bias.  In  this  document,  the 
parameter  bias  has  been  related  to  the  transplanting  of  parameters 
from  the  approximate  model  to  the  true  model.  Past  workers  have 
applied  an  ad  hoc  technique  that  partially  alleviates  the  noise  sen¬ 
sitivity  of  Prony's  method.  This  technique  consists  of  setting 
11  =  2n  (no  redundant  data)  and  setting  the  model  order,  n,  de- 


—  6(5— 


liberately  high  to  make  M  large.  Since  >1  is  the  number  of  samples 
used  in  the  estimation  procedure,  a  large  value  of  M  corresponds  to 
using  a  lot  of  data.  Hence,  this  ad  hoc  technique  uses  a  large 
amount  of  data  just  as  least-squares  Prony's  method  does  but  does 
not  suffer  from  the  parameter  bias  of  the  least-squares  method.  Un¬ 
biased  parameters  result  because  when  M  =  2n,  the  least-squares  pro¬ 
cedure  reduces  to  a  curve-fitting  procedure  which  must  be  unbiased. 

One  side  effect  of  this  technique  is  that  extra  "curve-fitting  poles" 
which  do  not  correspond  to  process  poles  are  introduced  in  the  desired 
exponential  representation . 

The  general  identification  scheme  presented  in  this  work 
opens  up  a  whole  range  of  possible  techniques  that  can  be  invented 
simply  by  making  different  choices  for  the  filters  within  the  identi¬ 
fication  model.  One  direction  in  which  future  research  may  be  fruit¬ 
ful  is  in  inventing  schemes,  other  than  least-squares  schemes,  to  com¬ 
bine  more  data  than  is  necessary  to  determine  the  model  parameters. 
Such  a  scheme  could,  perhaps,  be  realized  with  certain  choices  for 
the  filters  in  the  identification  model. 


APPENDIX 


Theorem:  If  G  is  a  singular,  complex  or  real,  NxN  dimensional 

matrix,  then  A.. A..  =  A..A..  for  ary  i  or  j ,  where  A  denotes  the 
13  31  11  jj  pq 

til  t^l 

element  of  the  p  row  and  the  q  ‘  column  of  adj  G  (that  is,  the 
cofactor  of  g  .) 

pq 


Proof :  Since  G  is  singular,  det  G=0.  Also,  from  the  identity 
I  det  G  =  G  adj  G,  it  follows  that  G  adj  G  =  0.  Each  column  of  adj  G 
then  must  be  a  solution  of  G_x  -  0-  Since  all  solutions  of  Gx  =  0 
are  collinear,  it  follows  that  the  columns  of  adj  G  are  also  collinear. 
Hence , 


adj  G  ** 


Vl  Xiy2 


Vl  V2 


*  •  • 


Vn 


Vn 


=  x  y 


(23) 


Vl  Vz  *  *  *  Vn 

T  T. 

where  x  =  [x^x9 — x^]  is  a  solution  of  Gx=0  and  y=  [y^y?...yjj]  is  the 

vector  required  for  equality  of  (23).  Then,  A^A^  =  (x^yj  (x^y^)  = 

(x.y.) (x.y.)  =  A. .A. .. 

11  3^3  11  3J 
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